VDEM <- read.csv("V-Dem-DS-CY-v6.1.csv")
VDEM <- subset(VDEM, select=c(COWcode, year, v2cltort, v2cltort_sd, v2clkill, v2clkill_sd, v2csreprss, v2csreprss_sd, country_name))
names(VDEM) <- c("COW", "YEAR", "v2cltort", "v2cltort_sd", "v2clkill", "v2clkill_sd", "v2csreprss", "v2csreprss_sd", "country_name")


# define objects to save model estimates for plotting
diffs <- NA
diffs.se <- NA
beta1 <- beta2 <- NA
se1 <- se2 <- NA


year <- 1949:2010


# loop through all the different model specifications for this treaty variable
for(k in 1:2){
for(j in 1:length(year)){
    
    data <- VDEM
    
    if(k==1)data$VAR <- data$v2cltort
    if(k==1)data$VAR.SD <-data$v2cltort_sd
    if(k==2)data$VAR <- data$v2clkill
    if(k==2)data$VAR.SD <-data$v2clkill_sd
    #data$VAR <- data$v2csreprss
    #data$VAR.SD <-data$v2csreprss_sd
    
    OBS <- rev(cumsum(rev(as.numeric(table(data$YEAR)))))

# load latent treaty data with latent treaty estimates included
treaty <- read.csv("CYEstimateLatentHRtreaty04.csv")

# create treaty proportion variable out of binary treaty variables
test <- treaty[,4:(ncol(treaty)-2)]
NA.test <- !is.na(test)

TREATY <- apply(test,1,sum,na.rm=T)
TOTAL <- apply(NA.test,1,sum,na.rm=T)


treaty$TREATY <- TREATY / TOTAL

treaty <- subset(treaty, select=c(ccode, year, TREATY))
names(treaty) <- c("COW", "YEAR","TREATY")

# load polity data
polity <- read.csv("p4v2010.csv", na.strings = c("-66", "-77", "-88", "-99"))
polity2 <- subset(polity, select=c(ccode, year, polity2))

# load updated Gleditsch trade gdp and population data
new.gdp <- read.delim("gdpv6.txt")
new.gdp <- subset(new.gdp, select=c(statenum, year, rgdppc, pop))
names(new.gdp) <- c("ccode", "year", "rgdpch", "pop") # rename because of the old name


# merge the data together and make lags
data <- merge(data, treaty, by.x=c("COW", "YEAR"), by.y=c("COW", "YEAR"), all.x=TRUE, all.y=FALSE, suffixes = c("",".treaty"))
nrow(data)


data <- merge(data, polity2, by.x=c("COW", "YEAR"), by.y=c("ccode", "year"), all.x=TRUE, all.y=FALSE)
nrow(data)

data <- merge(data, new.gdp, by.x=c("COW", "YEAR"), by.y=c("ccode", "year"), all.x=TRUE, all.y=FALSE)
nrow(data)


data.lag <- data
data.lag$YEAR <- data.lag$YEAR + 1
data <- merge(data, data.lag, by.x=c("COW", "YEAR"), by.y=c("COW", "YEAR"), all.x=TRUE, all.y=FALSE, suffixes = c("",".lag"))
#data$time <- data$YEAR-STARTYEAR
#data <- subset(data, !is.na(latentmean.lag) & !is.na(latentmean.treaty.lag) & YEAR>=1949 & YEAR<=2010)
#data <- subset(data, !is.na(latentmean.lag) & !is.na(latentmean.treaty.lag) & YEAR>=1949 & YEAR<=2005)
data <- subset(data, !is.na(VAR.lag) & !is.na(TREATY.lag) & !is.na(polity2.lag) & !is.na(rgdpch.lag) & !is.na(pop.lag) & YEAR>=year[j] & YEAR<=2010)
nrow(data)

#OBS <- rev(cumsum(rev(as.numeric(table(data$YEAR)))))

newdata <- list()

# take n draws from the posterior distribution and make n new datasets in the list object just created
for(i in 1:50){
  newdata[[i]] <- data
  newdata[[i]]$draw.lag <- rnorm(cbind(rep(1,nrow(data))), mean=data$VAR.lag, sd=data$VAR.SD.lag)
  newdata[[i]]$draw <- rnorm(cbind(rep(1,nrow(data))), mean=data$VAR, sd=data$VAR.SD)
}


# define model formula for each specification
FML <- "VAR ~ draw.lag + TREATY.lag"
#FML <- "VAR ~ draw.lag + TREATY.lag + polity2.lag"
#FML <- "VAR ~ draw.lag + TREATY.lag + polity2.lag + log(rgdpch.lag)"
#FML <- "VAR ~ draw.lag + TREATY.lag + polity2.lag + log(rgdpch.lag) + log(pop.lag)"
#FML <- "VAR ~ draw.lag + TREATY.lag + log(rgdpch.lag) + log(pop.lag)"
#FML <- "VAR ~ draw.lag + TREATY.lag + log(rgdpch.lag)"
#FML <- "VAR ~ draw.lag + TREATY.lag + log(pop.lag)"
#FML <- "VAR ~ draw.lag + TREATY.lag + polity2.lag + log(pop.lag)"

# define function to combine multe model estimates (this incorporates uncertatiny from the latent variables, see the supplementary appendix section F for more details about this function)
milm <- function(fml, midata){
 xx <- terms(as.formula(fml))
 lms <- matrix(data=NA, nrow=(length(attr(xx, "term.labels")) + 1), ncol=length(midata))
 ses <- matrix(data=NA, nrow=(length(attr(xx, "term.labels")) + 1), ncol=length(midata))
 vcovs <- list()
  for(i in 1:length(midata)){
    tmp <- lm(formula=as.formula(fml), data=midata[[i]])
    lms[,i] <- tmp$coefficients
    ses[,i] <- sqrt(diag(vcov(tmp)))
    vcovs[[i]] <- vcov(tmp)
  }
 par.est <- apply(lms, 1, mean)
 se.within <- apply(ses, 1, mean)
 se.between <- apply(lms, 1, var)
 se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(midata))))
 list("terms"=names(tmp$coefficients), "beta" = par.est, "SE"=se.est, "vcovs"=vcovs,"coefs" = lms)    
}

# call the milm function with the list of datasets newdata
model <- milm(fml=FML, midata=newdata)

# create object for printing during the loop and then print the results to screen
temp <- data.frame(model$terms, model$beta, model$SE, model$beta/model$SE)

print(paste("Model Output"))
print(temp)

print(as.numeric(temp[temp$model.terms=="TREATY.lag",2]))
print(as.numeric(temp[temp$model.terms=="TREATY.lag",3]))



#par(mar=c(5,7,5,5))

VDEM.est <- as.numeric(temp[temp$model.terms=="TREATY.lag",2])
VDEM.se <- as.numeric(temp[temp$model.terms=="TREATY.lag",3])


beta1[j] <- VDEM.est
se1[j] <- VDEM.se

print(c(VDEM.est, VDEM.se))
}






# barplot of the differences across all 8 model specifications
if(k==1) par(mfrow=c(2,1))
if(k==1) par(mar=c(2,4,4,6))
if(k==2) par(mar=c(4,4,2,6))
barplot(beta1, space=0, ylim=c(-.1,1.1), col=grey(.9), yaxt="n", border="white", xlim =c(1,60))
abline(h=0, col=grey(.5))
for(i in 1:7){abline(v=(c(2,12,22,32,42,52,62)-.5)[i], lwd=.75, col=grey(.75), lty=2)}
for(i in 1:length(year)){
    lines(c(i-.5,i-.5), c(beta1[i]+1.96*se1[i], beta1[i]-1.96*se1[i]), lwd=1)
    lines(c(i-.5,i-.5), c(beta1[i]+se1[i], beta1[i]-se1[i]), lwd=2)
}
#box()
#axis(side=1, at=c(.5:61.5), labels=rep("", 62))
#axis(2, at=c(-.20, -.15, -.10, -.05, 0, .05, .10, .15, .20), las=2)
if(k==1)mtext(side=4, "Model Coefficient \nV-DEM Torture DV", line=3.0, cex=1.0)
if(k==2)mtext(side=4, "Model Coefficient \nV-DEM Killing DV", line=3.0, cex=1.0)
if(k==1)mtext(side=3, "Human Rights Proportion \nTreaty Variable Coeffients and Differences in Coefficients", line=-1.5, cex=1.25)

axis(side=1, at=c(.5:61.5), labels=rep("", 62))
if(k==2) axis(1, at=c(2,12,22,32,42,52,62)-.5, label=c(1950, 1960, 1970, 1980, 1990, 2000, 2010), tick=F, cex.axis=1.5)
#axis(2, at=c(-.35,-.30, -.25, -.20, -.15, -.10, -.05, 0, .05, .10, .15, .20, .25, .30, .35), las=2)
axis(2, at=c(-.10, 0, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1, 1.1), las=2)

}

